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Q_; Abstract 

We provide a method for analytically constructing high-accuracy templates for the gravitational 

O _ 

. wave signals emitted by compact binaries moving in inspiralling eccentric orbits. By contrast to the 



simpler problem of modeling the gravitational wave signals emitted by inspiralling circular orbits, 
which contain only two different time scales, namely those associated with the orbital motion and 



. the radlatlou reaction, the case of inspiralling eccentric orbits involves three different time scales: 

o : 

■rj" , orbital period, periastron precession and radiation-reaction time scales. By using an improved 



'method of variation of constants', we show how to combine these three time scales, without 
making the usual approximation of treating the radiative time scale as an adiabatic process. We 
explicitly implement our method at the 2.5PN post-Newtonian accuracy. Our final results can 
be viewed as computing new 'post-adiabatic' short period contributions to the orbital phasing, 



o 
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X 

■ or equivalently, new short-period contributions to the gravitational wave polarizations, /i+.xj that 

should be explicitly added to the 'post-Newtonian' expansion for /i+.x , if one treats radiative effects 
on the orbital phasing of the latter in the usual adiabatic approximation. Our results should be of 
importance both for the LIGO/VIRGO/GEO network of ground based interferometric gravitational 
wave detectors (especially if Kozai oscillations turn out to be significant in globular cluster triplets), 
and for the future space-based interferometer LISA. 

PACS numbers: 04.30Db, 04.25.Nx, 04.80.Nn, 95.55.Ym 
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I. INTRODUCTION 



Inspiralling black hole binaries are considered to be the most probable source of detectable 
gravitational radiation for the first generation laser interferometric gravitational-wave de- 
tectors that are operational or nearing the completion of their construction phase [1]. The 
above understanding is based on both astrophysical considerations [2] and the availability of 
highly accurate general relativistic theoretical waveforms required to pluck the weak gravita- 
tional wave signal from the noisy interferometric data [3]. Inspiralling compact binaries are 
usually modeled as point particles in quasi- circular oihits. For long lived compact binaries, 
the quasi-circular approximation is quite appropriate, as the radiation reaction decreases 
the orbital eccentricity to negligible values by the epoch the emitted gravitational radiation 
enters the sensitive bandwidth of the interferometers. It is easy to deduce that for an iso- 
lated binary, the eccentricity goes down roughly by a factor of three, when its semi-major 
axis is halved [4]. 

In recent times, however, scenarios involving compact eccentric binaries are being sug- 
gested as potential gravitational wave sources even for the terrestrial gravitational-wave 
detectors. For instance, one such proposed astrophysical scenario [5], involves hierarchical 
triplets (say, 123), usually modeled to consist of an inner (say, 12) and an outer binary (say, 
03, where denotes the center of mass of the 12 binary). If the mutual inclination angle 
between the orbital planes of the inner and of the outer binary is large enough, then the time 
averaged tidal force on the inner binary may induce oscillations in its eccentricity, known in 
the literature as the Kozai mechanism [6]. It was shown that in globular clusters, the inner 
binaries of hierarchical triplets undergoing Kozai oscillations can merge under gravitational 
radiation reaction [5]. Later, it was shown that a good fraction of such systems will have 
eccentricity ~ 0.1, when emitted gravitational radiation from these binaries passes through 
10 Hz [7]. 

It is also believed that the effect of orbital eccentricity will have to be modeled accurately 
while computing theoretical waveforms for compact binaries relevant for the space-based 
laser interferometric gravitational- wave detector LISA [8] . The above statement is supported 
by a recent finding that by observing stellar mass black hole binaries in highly eccentric orbits 
- which may be common in globular clusters - one can estimate accurately not just the 
intrinsic binary parameters like masses and eccentricity, but even the position of the host 
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cluster [9]. It is also well known that the cosmological supermassive black hole binaries, 
embedded in surrounding stellar populations, would be powerful gravitational wave sources 
for detectors like LISA [10]. However, the question whether these binaries will be in a quasi- 
circular or an inspiralling eccentric orbit by the time their gravitational waves are detectable 
by LISA is not yet settled [11]. Recently, it was suggested that if the Kozai mechanism were 
operative, these supermassive black hole binaries, in highly eccentric orbits, would merge 
within the Hubble time [12]. Very recently, it was shown that, using very-long-baseline 
interferometer (VLBI), the unresolved core of the radio galaxy 3C 66B executes well-defined 
elliptical motions with an yearly period, which was interpreted as the first direct evidence 
for the detection of a supermassive black hole binary [13]. The above observation raises the 
interesting possibility of being able to detect gravitational waves from supermassive black 
hole binaries in eccentric orbits using LISA. 

Even though various versions of 'ready to use' high accuracy search templates for inspi- 
ralling compact binaries with arbitrary masses in quasi- circular oihit exist [3], so far none is 
available for compact objects in inspiralling eccentric orbits. Before characterizing the strat- 
egy and new results presented in this paper, let us summarize the relevant existing literature 
on the influence of eccentricity on the gravitational wave polarizations, and /ix, with and 
without including the gravitational radiation reaction. After the seminal work of Peters and 
Mathews [14], it was in the context of spacecraft Doppler detection of gravitational waves 
from compact binaries that the first explicit expressions for Newtonian accurate gravita- 
tional wave polarization states were derived [15]. Using the method of osculating orbital 
elements and a numerical integration approach the effects of eccentricity and dominant radi- 
ation damping on and was studied in [16]. First- and First-and-half-post-Newtonian 
accurate analytic expressions for far-zone fluxes and gravitational wave polarizations, for 
compact binaries in eccentric orbits, were computed in a series of papers in [17]. Using 
Newtonian accurate orbital motion, Refs. [18, 19] studied the effect on gravitational wave 
polarizations of introducing by hand some secular effects either in the longitude of the perias- 
tron [18] or in the semi-major axis and eccentricity [19]. Using such waveforms, the influence 
of eccentricity on the signal to noise ratio in gravitational wave data analysis, was examined 
in [19-21]. These waveforms were also used to show that LISA will be sensitive to eccentric 
Galactic binary neutron stars [22] and that by measuring their periastron advance, accurate 
estimates for the total mass of these binaries may be obtained [23]. However, the widely 
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used gravitational wave templates, to detect gravitational waves from compact binaries in 
quasi- circular orbits, are based on the second post-Newtonian accurate expressions for h+ and 
/ix, supplemented by expressions giving adiabatic time evolution for the orbital phase and 
frequency also to the second post-Newtonian order [24, 25] [However, see [26] for the (numer- 
ical) construction of gravitational wave templates going beyond the adiabatic approximation 
in the case of quasi- circular orbits]. The second post-Newtonian order, usually referred to 
as the 2PN order, gives corrections to leading order contributions in gravitational theory, to 
order of ~ {v/c)'^ ~ {Gm /c^r)^, where m,v and r are the total mass, orbital velocity and 
the separation of the binary. The 2PN contribution to the gravitational waveform, required 
for the construction of and h^, and the associated far- zone fluxes for binaries moving in 
general eccentric orbits, in harmonic gauge, were computed in [27, 28]. Employing the 2PN 
accurate generalized quasi-Keplerian parametrization of Damour, Schafer and Wex [29, 30] 
available in Arnowitt, Deser and Misner (ADM) coordinates, to represent relativistic motion 
binaries in eccentric orbits, 2PN corrections to the rate of decay of the orbital elements of the 
representation as well as the explicit expressions for /i+ and hx were provided in [28, 31, 32]. 
The above mentioned expressions for /i+ and hx represent gravitational radiation from an 
eccentric binary, during that stage of inspiral where the gravitational radiation reaction is so 
small that the orbital parameters can be treated as essentially unchanging over a few orbital 
periods ('adiabatic approximation'). The effects of eccentricity, advance of periastron and 
orbital inclination on power spectrum of the dominant Newtonian part of the polarizations 
were also presented in [31]. 

The aim of this paper is to provide a method for explicitly constructing high-accuracy 
waveforms emitted by compact binaries moving in inspiralling eccentric orbits. Compared 
to the existing high-accuracy waveforms for inspiralling circular orbits [3], the inclusion 
of orbital eccentricity into such templates is a non-trivial task as eccentricity brings along 
a new physical aspect, the precession of the periastron, and thus one must contend with 
precession and radiation reaction at the same time. In the quasi-circular case, there are two 
time-scales related to the orbital period and the radiation reaction. In the quasi-eccentric 
case, one has a third additional time-scale related to the precession of periastron. The 
technical problem we shall tackle and solve below is that of combining, in a consistent 
framework these three time scales, without making the usual approximation of treating 
the radiation-reaction time scale as an adiabatic process. We then explicitly implement 
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our method at the '2PN+2.5PN' accuracy, i.e. the effect of perturbing a '2PN-accurate 
analytic' description of eccentric orbits by the 2.5PN level radiation-reaction. It is useful to 
note that the gravitational wave observations of inspiralling compact binaries, are analogous 
to the high precision radio-wave observations of binary pulsars. The latter makes use of an 
accurate relativistic 'timing formula' which requires the solution to the relativistic equation 
of motion for a compact binary moving in an elliptical orbit while the former demands 
accurate 'phasing', i.e. an accurate mathematical modeling of the continuous time evolution 
of the gravitational waveform. The mathematical formulation, which resulted in an accurate 
'timing formula', as given in [34, 35], was obtained by one of us many years ago [36, 37]. 
The present work will rely on techniques from [36, 37] to combine the mentioned above three 
time scales to implement post-Newtonian accurate 'phasing' for compact binaries moving in 
elliptical orbits. 

The plan of the paper is the following. In Sec. II, we outline the steps required to do the 
'phasing'. The method we use to find the domain of validity of our approach is presented in 
Sec. III. In Sec. IV, we formulate the procedure required to construct time evolving hx and 
We apply, in Sec. V, the formalism to do the 2.5PN accurate phasing. Sec. VI displays 
the expressions required to extend the 'phasing' to higher PN orders. Finally, in Sec. VII, 
we summarize our approach, results and point out possible future extensions. 

II. THE PHASING OF GRAVITATIONAL WAVEFORMS 

The theoretical templates for compact binaries, required to analyze the noisy data from 
the detectors, as noted earlier, usually consist of /i+ and h^, the two independent gravita- 
tional wave (GW) polarization states, expressed in terms of the binary's intrinsic dynamical 
variables and location. These two basic polarization states /i+ and are given by 



where hjj^, the transverse-traceless (TT) part of the radiation field, is expressible in terms 
of a post-Newtonian expansion in (v/c), and where p and q are two orthogonal unit vectors 
in the plane of the sky i.e. in the plane transverse to the radial direction linking the source 
to the observer. 
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The TT radiation field is given, by the existing gravitational wave generation formalisms 
[24, 25, 27], as a post-Newtonian expansion of the form 



hi + -hi + \hl + + \hi + \hl + \h% + --- 



where, for instance, the leading ('quadrupolar') approximation is given (in a suitably defined 
'center of mass frame', see below) in terms of the relative separation vector x and relative 
velocity vector v as 

Gm 
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where Vijkmi^) is the usual transverse traceless projection operator projecting normal to 
N, where N = R'/i?', R' being the radial distance to the binary. The reduced mass of 
the binary fi is given by mim2/m, where m = rrii + m2 is the total mass of the binary 
consisting of individual masses mi and m2. We also used Vij = ViVj, and Uij = riirij, 
where Vi and are the components of the velocity vector v = d'^/dt and the unit relative 
separation vector n = x/r, where r = |x|. When inserting the explicit expression of and 
its higher-PN analogues h\j, hfj ■ ■ ■ which are currently known up to hfj [27, 28, 31], {h^j is 
currently being explicitly derived [38]), one ends up with a corresponding expression for the 
two independent polarization amplitudes, Eqs. (1), as functions of the relative separation r 
and the 'true anomaly' 0, i. e. the polar angle of x, and their time derivatives. 
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For instance, if we follow the conventions used both in [27] and [25], for choosing the or- 
thonormal triad p, q, N (namely, N from the source to the observer and p toward the 
correspondingly defined 'ascending' node), i.e. if we use 



x = p r cos + (q cos i + N sin i)r sin i 



(5) 



where i denotes the inclination of the orbital plane with respect to the plane of the sky, the 
lowest-order contribution in Eq. (4) reads 



\hl{r,(l),r, 
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c4 ^^ '^' '^^ c'^R 
where rj = fi/m = mim2/ {nii + 7712)^, C and S are shorthand notations for cos 2 and 
sinz. The orbital phase is denoted by 0, = d(f)/dt and r = dr/dt = n ■ v, where v = 
p (r cos (j)—r(j) sin 0) + (q cosi+N sini) (r sin 0+r cos 0). We note that in our expressions 
for /ix and the coefficients of cos 2 and sin 2 differ from those derived in [31] by a 
minus sign as in that paper the true anomaly was measured from q, rather than from the 
line of ascending node as done here and in [25]. 

Having in mind the existence of expressions such as Eq. (6) giving the GW amplitudes 
hx in terms of the relative motion x, v of the binary, it is clear that they must be 
supplemented by explicit expressions describing the temporal evolution of the relative mo- 
tion, i.e. describing the explicit time dependences r(t), 0(t), r(t), and 0(t). We refer to 
as phasing, any explicit way to define the latter time-dependences, because it is the crucial 
input needed beyond the 'amplitude' expansions, given by Eqs. (4), to derive some ready to 
use waveforms h+^xit). 

Let us note two things about the structure sketched above for /i+,x- First, the possibility 
to express the GW polarization amplitudes in terms of only the relative motion x, v, relies 
on the possibility to go to a suitable center-of-mass frame. The validity of a centre-of-mass 
theorem up to order inclusive i. e. in the presence of the leading radiation reaction was 
first shown in [36] (in harmonic coordinates). The analogous result, in ADM coordinates, 
was obtained in [39], where it was shown that there existed six first integrals of the 2.5PN 
equations of motion: a total momentum "P^^g-j = "^(0) + ^"^"^(2) +'^"^'^(5) 'center- 

of-mass' constant /C*^^^ = Q^^^-^ — ^"^(5) , which could both be set to zero by applying a suitable 
Poincare transformation. The recent obtention of (manifestly or not) Poincare invariant 3PN 
equations of motion [40, 41] and the construction of a corresponding complete set of 3PN 
conserved quantities [42, 43] allows one to extend the construction of and /C* = — tP* 
to order inclusive and thereby define a 3PN-accurate center-of-mass frame [44]. [Note, 
however, that the contribution ^^^^^ to Q^^q^ introduced by [44] coincides with ^^^^^ of 
[36] only in the center-of-mass frame.] At the next PN level, the above 3PN "conserved 
quantities" will not be conserved anymore because of the 3.5PN component of radiation 
reaction so that 
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^(6) = ^(^) ' Re) = 0{^) . (7) 
The (9(c~^) 'recoil' of the center-of-mass imphed by Eqs. (7) is expected to influence the 
waveform only at the 0{c~^) level. Indeed, if we think of the binary as a GW source emitting 
the 'relative' signal, given by Eqs. (6), in its (instantaneous) rest frame (namely, the above 
defined 3PN center of mass frame), the time-dependent recoil of the latter rest frame will 
introduce both a N ■ vcm/c Doppler shift of the phasing and a corresponding modification 
of the amplitudes /?.+,x- 

Second, we should mention that the possibility to express /i+^x only in terms of r, 
and their time derivatives holds because we restricted ourselves to non-spinning objects. 
In the presence of spin interactions, the orbital plane is no longer fixed in space and one 
needs to introduce further variables, notably a (slowly varying) 'longitude of the node' Q. 
Correspondingly, the polarization direction p cannot be defined anymore as the line of nodes. 
We note in this respect that such a situation was dealt with in the problem of the timing 
of binary pulsars [35] and it might be advantageous to use the conventions used there to 
define p and q. Namely, in terms of Fig. 1 of [35], p = Iq, q = Jo, but note that the 
binary pulsar convention uses as the third vector Iq x Jq, the direction from the observer to 
the source. Such a convention being natural when one thinks of the actual observation of a 
signal somewhere on the sky, as seen by us! 

Finally, one should make clear the coordinate systems that we shall use. Indeed, the 
explicit functional forms for h_^.{r,(f),r,(f)),hx{r,(j),r,(j)) as well as the phasing relations 
r{t),(f){t),r{t) and (j)(t) depend on the coordinate system used, though the final results 
hj^{t) and hx{t) do not [Note that /i^-^ and therefore h+{t) and hx{t) are coordinate in- 
dependent asymptotic quantities ]. Here, one has to face a slight mismatch between the 
harmonic coordinate systems in which standard GW generation formalisms derive the am- 
plitude expressions Eqs. (4) above, and the ADM coordinate systems which allow one to 
derive (when neglecting radiation reaction) a rather simple and elegant 'quasi-Keplerian' 
form of the general (eccentric) orbital motion [29, 30], and therefore of the phasing of the 
GW signal. As our work is focused on the latter phasing issue (in presence of radiation reac- 
tion), we shall consistently work in ADM- type coordinate systems because they will allow us 
to write down explicit analytical expressions for the orbital phasing r{t) and 0(t). We shall 
then assume that the starting amplitude expressions are first transformed from the original 
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harmonic-coordinates form /i+^x (yi? y2, yi? y2) to the corresponding ADM-coordinates ones 
/i+^x (xi, X2, xi, X2). Note that we do mean expressing /i+^x in terms of ADM positions and 
velocities, though the use of positions and momenta a priori looks more natural in the ADM 
Hamiltonian framework. All the formulas necessary for the transformation between the two 
systems have been worked out, at the 2PN order in [29, 45] and in [42, 43] for the 3PN. 
Note that the reduction to the center-of-mass can equivalently be performed before or after 
the transformation (ya, ya) (xa,Xa) [44]. Evidently, this transformation, which starts at 
2PN order, does not affect the lowest-order expressions exhibited above, Eq. (6). 

We have explicitly computed, in ADM coordinates, h^ ^^, ^+,x) ^+,x &nd /i+^x; which give 
FN corrections to /i+.x to 2PN order, in terms of r, 0, f, 0, in the convention used to obtain 

X . In appendix A, we briefly describe the steps to get the above corrections and sketch 
the structural forms of these corrections. 

III. DELINEATING THE 'STABLE' ECCENTRIC ORBITS AND THE 'QUASI- 
KEPLERIAN' ONES 

In the case of inspiralling czrcw/ar orbits, a very important role, for the phasing of gravita- 
tional waves, is played by the last stable (circular) orbit (LSO). In the zeroth approximation, 
circular orbits above the LSO, and in the presence of radiation reaction, can be described as 
an adiabatic sequence of circular orbits. This approximation breaks down when the binary 
reaches the LSO, at which point the orbital motion changes into a kind of (relative) plunge. 
To describe the smooth transition between the adiabatic inspiral and the plunge, one needs 
to use a formalism for the orbital motion (such as the 'effective one body' approach [46]), 
which goes beyond the usual, purely perturbative, post-Newtonian approach. 

In the present paper, we study inspiralling eccentric orbits, evolving under radiation 
reaction, and we shall use a purely perturbative post-Newtonian approach (but one which 
goes beyond the zeroth order, adiabatic approximation). Such a treatment can only be valid 
if we stay sufficiently above any eccentric analog of the LSO, i.e. if we consider eccentric 
orbits which are 'stable' in the sense of being separated by a potential barrier from any 
plunge motion. The purpose of this section is to delineate, in the plane of the parameters 
representing the two-dimensional manifold of eccentric orbits, the locus where such orbits 
would cease to be bona fide bound orbits to become plunge-type ones. To do this, we need 
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to use the non-perturbative effective one body (EOB) formahsm, because the transition 
between bound and plunge orbits has a non-perturbative, strong-field origin. As the rest 
of the paper will consider the conservative part of the orbital motion at the second post- 
Newtonian (2PN) approximation, we shall consistently use the EOB Hamiltonian at the 
resummed 2PN level [46]. (To work at the 3PN level one should use the more complicated 
EOB formalism derived in [47]). 

The real (resummed 2PN) EOB Hamiltonian describing the conservative part of the 
orbital motion {i.e. when neglecting radiation reaction), has the form (we recall that rj = 
fj,/m = mim^l {mi + m2Y) 



Hreal{R, Pr, P<t>) 



mc 



'eff 



1 ) , where. 



cff 



cff 



'A{R) 
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PI 
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A{R) 



2Gm 



2r] 



Gm 



(8) 



(9) 



(10) 



c'^R ' \d^R^ 

Here R denotes the effective radial separation between the two bodies, Pr the corresponding 
(relative) radial momentum, and J7 the (relative) total angular momentum. The total 
energy (including the rest mass) will be denoted by Sj-ca\ (or simply S when no confusion 
with the 'effective' energy can arise ). The total energy Sj-eai = -f^reai is related to the 

^2 'l + 2r/(Ee 



'effective' specific energy E^fi = H^s by Sj-c^i = mc \jl + 21] ( E^fi — 1 ). We shall not need 
the explicit expression of the EOB metric component B{R), but only use the fact that 
B{R) > 0. Taking into account the fact that the radial kinetic energy term Pl/{^i^c^B{R)) 
in Eq. (9) is positive, the radial motion can be qualitatively understood in terms of the 
angular-momentum dependent (effective) 'radial potential' 



Wj{R) = A{R) 



1 + 



(11) 



which is a generalization (to the comparable mass case) of the well-known radial potential 
for test-particle orbits around a Schwarzschild black hole, namely. 



W^{R) = (1 - 2Gm/c'R){l + j' /{fi'c'R')) 



(12) 



(which is simply the ?7 — limit of Eq. (11)). 
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In Fig. 1, we present typical plots of Wj{R) on which we mark 'energy levels' correspond- 
ing to some eccentric and circular orbits. This plot makes it clear that an a priori bound 
motion {i.e. with total energy £^rcai < rnc^, i.e. E^s < 1) will execute some precessing but 
stable motion only if the energy- level H^f^ = E^^ stays above the radial potential Wj{R) in a 
finite radial interval [-Rmin 7 -Rmax]- The 'centrifugal barrier' preventing this radial motion to 
plunge towards smaller separations is on the left of -Rmin- Therefore the locus of orbits which 
are on the verge of plunging corresponds to the case where the energy-level (horizontal) line 
would be tangent to the top of the potential barrier, i.e. to 



In the domain of interest, the Eq. (13), considered for some given JT, will have two roots, 
say Rp{J) and Rc{J), with Rp < Re- The larger root Rc{J) defines the set of stable circular 
orbits. Note, in this respect, that the smaller root Rp{J), which marks the 'plunge' locus 
also corresponds to the locus of unstable circular orbits. Finally, the domain in the energy- 
angular momentum plane corresponding to bound (non plunge) motion is restricted by the 
inequalities 



which defines, when using the link (8), a corresponding double inequality involving J and 
£^rcai- Note that, as one approaches the plunge boundary, i.e. for orbits close to the orbit 
marked 'p' in Fig. 1, the character of the orbital motion starts to deviate very much from that 
of a usual, perturbative, slowly precessing, "quasi-Keplerian" motion. Instead, it becomes 
what is referred to as a "zoom- whirl" motion in [48], i.e. a motion which alternates between 
one large-excursion elliptic-Keplerian-like orbit and several quasi-circular orbits near the 
periastron (which, as noted above is close to an unstable circular orbit). As the formalism 
we shall use in this paper to analytically represent the orbital motion assumes a "quasi- 
Keplerian" representation (see below), we shall need to stay sufficiently away from the 
plunge boundary to ensure the numerical validity of such a representation. Before coming 
to the issue of what one exactly means by "sufficiently away", let us finish describing the 
analytical estimate of the plunge boundary, as defined by the inequalities given in Eq. (14). 

Let us analytically estimate the two crucial roots Rp{J) and Rc{J) of Eq. (13). Using 
the dimensionless, scaled variables u = Gm/c^R, j = cj /{jiGm), the radial potential reads 



dWj{R) 
dR 



0. 



(13) 



Wj{R,{J)) < El^ < Wj{Rp{J)) , 



(14) 
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W{u) = A{u) (1 +jV) (15a) 
where, A{u) = 1 - 2u + 2r]u^ . (15b) 



3u^ -u + — -r]u'^{— + 5u^) = 0. (16) 



The equation (13), or better —W'{u)/{2j^), reads 

1 2.3 
3 3 

When r] 0, Eq. (16) becomes a quadratic equation, with the two roots 

1 



4(J) g 



(17) 



where the plus sign corresponds to the 'plunge' boundary (larger u, i.e. smaller R), while 
the minus sign corresponds to circular orbits. An accurate (at least when > 12) analytical 
estimate of the 77- deformations of the above two roots, i.e. the roots of the quartic equation, 
Eq. (16), corresponding to Rp and Rc, is obtained by inserting expressions (17) into the 
ry-dependent terms in Eq. (16). This yields 



l±yi-^[l-^(4P(3 + 5j2C4)2)] 



(18) 



We have verified that the analytical estimate, Eq. (18), is a numerically accurate estimate of 
the two roots Up{j), Uc{j). Inserting this result (with = Up, u~ = Uc) into Eq. (14) then 
yields an explicit (2PN-level) estimate of the domain of 'non-plunge' eccentric orbits in the 
{S, J') plane. Another way to describe the plunge boundary in the {S, J') plane, which does 
not need to assume that > 12, is to give a parametric representation of this boundary 
in terms of the parameter u in the form S = S{u),J' = J{u). This is simply obtained by 
solving Eq. (16) for j, which gives j{u) = y/^-^yu ^ then obtains £ = £{u) by 

W u (1— 3n+5 r] u'^) 

substituting for J by in Yiqs. (8) and (9). 

Having discussed the location of the plunge boundary, let us now discuss the issue of 'how 
far from the boundary' we need to stay for allowing us to rightfully use the analytical quasi- 
Keplerian representation of [29, 30] (see Eqs. (26) below). As this representation assumes, 
among other things, that the orbits are slowly precessing, we need to set an upper limit on the 
rate of periastron precession. [This will then, de facto, eliminate the possibility of whirl- zoom 
orbits, which contain, during part of their orbital period, a rapidly precessing quasi-circular 
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motion]. To see more precisely what upper limit we should set on periastron precession, 
let us go back to the EOB representation of the motion. From previous work on the 2PN- 
accurate EOB dynamics [46], we know that the comparable mass case (77 ~ 0.25) is rather 
close to the test-mass one [rj ^ 0), yielding geodesic motion in a Schwarzschild spacetime. 
Therefore, let us consider the domain of parameter space for which periastron precession 
around a Schwarzschild black hole is well described by a slowly precessing, quasi-Keplerian 
motion. For this, we consider the exact formula, given by Eq. (A8) in the Appendix A of [29], 
which gives the angle of return to the periastron for a test particle moving in Schwarzschild 
spacetime. It can be easily checked that, for the elliptical orbits ( the eccentricity parameter 
< 1) we are interested in, the term whose expansion is most slowly convergent is the 
prefactor (1 — p)"^/"^ in Eq. (A8) of [29]. We must therefore impose i| <^ 1 to have a 
slowly precessing, quasi-Keplerian motion. When this inequality is satisfied (together with 
< < 1), we expect that the 2PN- accurate expressions for n = T being the radial 
orbital period, and ef in terms of S and J', as derived in [29, 30], to be numerically accurate. 
In terms of dimensionless non-relativistic energy per unit reduced mass E = (S — m c^) / n 
and j, defined earlier as cJ'/{fiGm), the expressions for n and read [49] 

^ = (-2 E)^/2|l -- (15 -r/) (-2 i?) + ^^^ (555 + 3077 + 11^/2) 
el = 1 + 2E f + EU{l-ri) + {17-7 r])Ef\ + \2{2 + r] + 57]^) 



-(17 - 11^)-^ + (112 - 477/ + IQrf) E^f 



(-2 E"|3/2 

-3{5-27]){l + 2Efy (19b) 



J 

Using the above expressions one can approximately express ^ in terms of ^ and et, and 
define 

12 f^/^ 
e=-~12-^^ (20) 
3^ (l-e*) 

We can now specify 'what small means' in terms of e to ensure a decent convergence of 
the crucial factor (1 — e)~^/^ entering the periastron precession expression. A minimal 
requirement would be to impose e < |- Indeed, when e = | the 2PN-accurate expression for 
the periastron-advance parameter k' = (1 — e)~^^^ — 1, (which yields the periastron advance 
of nearly circular orbits) namely k' = j + ^ gives the exact value to an accuracy ~ 3%. 
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Choosing such a threshold, e < |, for 'staying sufficiently away from the plunge boundary' 
leads to the following constraint on the parameters ^ and ef. 



(l-l?)3/2 (12) 




Later, when we evolve orbital elements and gravitational waveforms, we make sure that 
the eccentric orbits we study lie inside this domain, defined by the above inequality. Let 

us emphasize that this restriction is due to our use of, in the next section, the generalized 
quasi- Keplerian representation, given by Eqs. (28), (29) and (30). We could go beyond the 
limit, given by Eg. (21), by using, instead of the generalized quasi- Keplerian representa- 
tion, the exact Schwarzschild-like motion (analytically expressible in terms of rather simple 
quadratures) in the FOB metric. This will be tackled in the near future. 

IV. A METHOD OF VARIATION OF CONSTANTS 

In this section, we introduce a version of the general Lagrange method of variation of 
arbitrary constants, which was employed to compute, within general relativity, the orbital 
evolution of the Hulse- Taylor binary pulsar [36, 37]. The method begins by splitting the 
relative acceleration of the compact binary A into two parts, an integrable leading part ^0 
and a perturbation part. A' as 



In this work, we will work at 2.5PN accuracy and accordingly choose ^0 to be the acceleration 
at 2PN order and A' to be the (leading) contribution to radiation reaction. It will, 
however, be clear that our method is general and can be applied, for instance, to a 3.5PN- 
accurate calculation where ^0 would be the conservative part of the 3PN dynamics, and 
A' the C(c~^) + (9(c~^) radiation reaction. The method ffist constructs the solution to the 
'unperturbed' system, defined by 



A = Ao + A' . 



(22) 



X = 



V 



(23a) 



V = A(x,v). 



(23b) 



The solution to the exact system 



X = 



(24a) 
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V = ^(x,v) , (24b) 

is then obtained by varying the constants in the generic solutions of the unperturbed system, 
given by Eqs. (23). The method assumes (as is true for Af^^"''^"'') that the 

unperturbed system admits sufficiently many integrals of motion to be integrable. For 
instance, if we work with Aq = A2PN, we have four first integrals: the 2PN accurate energy 
and 2PN accurate angular momentum of the binary. We denote these quantities, written in 
the 2PN accurate center of mass frame as ci and d^'- 

Ci = ^(xi,X2, Vi, V2)|2PN CM , (25a) 
4 = j7'i(xi,X2, Vi, V2)|2PN CM, (25b) 

with corresponding 3PN definitions of ci and c^, if we were working with = ^s'p^'^™*'™. 

The vectorial structure of c^, indicates that the unperturbed motion takes place in a 
plane. The problem is restricted to a plane even in the presence of radiation reaction [36]. 
We may therefore introduce polar coordinates in the plane of the orbit r and such that 
X = ircos0 + jrsin0 with, say, i = p, j = qcosi + Nsini (see above). The functional 
form for the solution to the unperturbed equations of motion, following [29, 36], may be 
expressed as 

dS 

r = S{l;ci,C2) ; f = n— (/; ci, Ca) , (26a) 

dW 

4, = X + W{l-CuC2) ; ^= {l + k)n + n-^{l;ci,C2), (26b) 

where and / are two basic angles, which are 27r periodic and C2 = \c2\- The functions 
S{1) and W{1) and therefore ^(/) and ^(0 periodic in / with a period of 2n. In the 
above equations, n denotes the unperturbed 'mean motion', given by n = ^, P being the 
radial (periastron to periastron) period, while k = A$/27r, A$ being the advance of the 
periastron in the time interval P. The explicit 2PN accurate expressions for P and k in 
terms of Ci and C2 are given in [29]. The corresponding 3PN accurate ones are given in [42]. 
The angles / and A satisfy, still for the unperturbed system, / = n and A = (1 + k)n, which 
integrate to 

I = n{t- to) + ci , (27a) 



^ We denote by A, the variable denoted by m in Ref. [36]. In most current literature including this one, m 
denotes the total mass of the binary. Also note that, the variable k here is related to the fc-variable in 
[31] say feci by /c = fcci/c^ 
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A = (1 + k)n{t - to) + cx , (27b) 

where to is some initial instant and the constants q and cx, the corresponding values for / 
and A. Finally, the unperturbed solution depends on four integration constants: Ci, C2, Q 
and Cx- 

At the 2PN order, one can write down explicit expressions for the functions S{1) and 
W{1). Indeed, the generalized quasi- Keplerian representation [29, 30] yields: 

S{l;ci,C2) = 0^,(1 — e,. cos n) , (28a) 
W{l;ci,C2) = {l + k)(v-l) + fysm2v + ^sm3v, (28b) 

where v and u are some 2PN accurate true and eccentric anomalies, which must, in Eq. (28), 
be expressed as functions of I, ci, and C2, say as v = V(/;ci,C2) = V{U{1] 01,02)) and 
u = U{1; 01,02). In the above equations, a,, and are some 2PN accurate semi-major 
axis and radial eccentricity, while and are certain functions, given in terms of ci and 
02- [To avoid introducing new notation, the eccentric anomaly is denoted by u following 
standard convention. It should not be confused with u = Gm/o^R employed in Sec. III. 
A similar comment applies to the function v below in the quasi-Keplerian representation 
and the magnitude of the relative velocity v] The definitions of 2PN accurate functions 
u = U{1; 01,02) and v = V{u) are available in [29, 30]. First, the function v = V{u) is 
defined by 

V = V{u) = 2 arctan (^^^^) f ) ' 

Second, the function u = IA{1) is defined by inverting the following 'Kepler equation' I = l{u) 

Z = n — sin u + ^ sin \^(n) + ^ iy{u) — u) . (30) 

Then the function v = V{1) is obtained by inserting u = U{1) in = V{u), i.e. V{1) = 
V{U{1)). Here Ct and are some time and angular eccentricity and ft and gt are certain 
functions of Ci and 02, appearing at the 2PN order. In our computations, we use the following 
exact relation for v — u, which is also periodic in u, given by 

V — u = 2 tan ^ , (31) 

V 1 — cos u I 



1- /l-e^ 

where (3s = — -. We note that the extension of such a generalized quasi-Keplerian 

representation to the 3PN order is easily possible when working in ADM type coordinates. 
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Let us now turn to use the explicit unperturbed solution, Eqs. (26) and (27), for the 
construction of the general solution of the perturbed system, Eqs. (24). This is done by 
keeping exactly the same functional form for r, r, and 0, as functions of / and A, Eqs. (26), 
i.e. by writing 

dS 

r = S{l;ci,C2) ; r = n—{l; a, C2) , (32a) 

dW 

(f) = X + W{l;ci,C2) ; (j) = (1 + k)n + n-^{l; ci, C2), (32b) 

but by allowing temporal variation in Ci = Ci(t) and C2 = C2(t) [with corresponding temporal 
variation in n = n{ci, C2) and k = k{ci, C2)], and, by modifying the unperturbed expressions, 
given by Eqs. (27), for the temporal variation of the basic angles / and A entering Eqs. (32) 
into the new expressions: 

I = / ndt + ci{t), (33a) 

J to 

ft 

A 

'to 



[ {l + k)ndt + cx{t) , (33b) 

Jto 



involving two new evolving quantities Q(t), and Cx(t). In other words, we seek for solutions 
of the exact system, Eqs. (24), in the form given by Eqs. (32) and (33) with four 'varying 
constants' Ci(t), C2(t), ci{t) and cx(t). The four variables {ci,C2,q,ca} replace the original 
four dynamical variables r, r, and 0. It can be verified that the alternate set {ci, C2, Q, cx} 
satisfies, like the original phase-space variables, first order evolution equations [36, 37]. These 
evolution equations have a rather simple functional form, namely, 

dr 

-^ = F^il-Cfs);a,P = l,2,l,X, (34) 

where the right hand side is linear in the perturbing acceleration. A'. Note the presence of 
the sole angle / (apart from the implicit time dependence of cp) on the RHS of Eqs. (34). 
The explicit expressions for these evolution equations were derived in [37], which in our 
notation read 

A'' , (35a) 



dci 


_ 9ci(x,v) 


~dt 




dC2 


_ 9C2(X,V) 






dci 




dt 





A" , (35b) 

^ / dS dci dS dc2 \ 
' (35c) 



dci dt dc2 dt 
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dcx dW dci dW dci dW dc2 
dt dl dt dci dt dc2 dt 

The evolution equations for ci and C2 clearly arise from the fact that Ci and C2 were defined 
as some first integrals in phase-space, say Eqs. (25). As shown in [37], there is an alternative 
expression for which reads 

dci f dQ\~^ f , dQ dci dQ dc2 



dt~[dl) V ■ 9ci dt dc2 dt ' ' ^^^^ 



where Q'^{1, Ci, C2) = r^(S'(/, Ci, C2), Ci, C2) and ^ is defined by 

dQ P dQ^ 



(37) 



dl An dr ' 

Both expressions, Eqs. (35c) and (36), involve formal delicate limits of the 0/0 form, for 
some (different) values of /. Taken together, they prove that these limits are well-defined and 
yield for dci/dt, an everywhere regular function of /. Anyway, the algebraic manipulation of 
the explicit forms of both Eqs. (35c) and (36) lead to well-defined expressions [For example, 
in the case of Eq. (35c), the problematic sinn factor in dS/dl, see Eq. (48a) below, nicely 
simplifies with a sin u factor present in the term within parenthesis on the RHS of Eq. (35c)] . 

The definition of / given by / = J^^ n{ca{t)) dt + Q(t), is equivalent to the differential 
form, ^ = + ^ = n + Fi{l,Ca):,a = 1,2, which allow us to define a set of differential 
equations for as functions of / similar to Eqs. (34) for Cq as functions of t. The exact form 
of the differential equations for Ca{l) reads 

dCa F(^{1^ Co) (r>Q\ 

~dr~n{ca) + Fi{i-cay ^ ' 

where Ca, a = 1,2, stands for ci and C2. Neglecting terms quadratic in F^, i.e. quadratic in 
the perturbation A! (e.g. neglecting (9(c^^°) terms in our application), we can simplify the 
system above to 

dcn 1 



- , ,F^{l;Ca) = Ga{l;Ca) ;a = l,2, 1, X; a = 1,2. (39) 
dl n[Ca) 

From here onward, we will neglect these 0{c~^^) terms in the evolution equations for €„(/), 
i.e. work with the simplified system, namely Eq. (39). At this stage, it is crucial to note, 
not only that the RHS of Eq. (39) is a function of Ci, C2 and the sole angle / (and not of 
A), but that it is a periodic function of /. This periodicity, together with the slow [Ga oc 
Fa (x A' = 0{c~^)] evolution of the c^'s, implies that the evolution of Ca{l) contains both 
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a 'slow' (radiation-reaction time-scale) secular drift and 'fast' (orbital time-scale) periodic 
oscillations. For the purpose of phasing, to model the combination of slow drift and the fast 
oscillations present in c^, we introduce a two-scale decomposition for Cq,(/) in the following 
manner 

Ca{l) = Ca{l) + Ca{l) , (40) 

where the first term Cq(/) represents a slow drift (which can ultimately lead to large changes 
in the 'constants' Ca) and Ca{l) represents fast oscillations (which will stay always small, 
i.e. of order 0{Ga) = 0{c~^)). This is proved by first decomposing the periodic functions 
Ga{l) (considered for fixed values of the other arguments ) into its average part and its 
oscillatory part: 

= ^ [ dlG{l,Ca), (41a) 

^TT Jo 

G^{l-Ca) = G^{l;Ca) -Ga{Ca). (41b) 

Note that, by definition, the oscillatory part Ga{l) is a periodic function with zero average 
over /. Then assuming that Ca in Eq. (40) is always small (c^ = 0{Ga) = 0{c~^)), one can 
expand the RHS of the exact evolution system, given by Eqs. (39), as 

^ + ^ = G„(/; + Ca) = G^il; Ca) + 0{Gl) , 

= G^il;ca) + G^{l;ca) + 0{GI). (42) 

We can then solve, modulo 0{G'^), the evolution equation (42) by defining Ca(/) as a solution 
of the 'averaged system' 

^ = G^ica) , (43) 
and by defining Cq,(/) as a solution of the 'oscillatory part' of the system 

^ = G^{l,Ca). (44) 

During one orbital period (0 < / < 27r) the quantities on the RHS of Eq. (44) change only 
by 0{Ga). Therefore, by neglecting again terms of order 0{G^) ~ (9(c~^°) in the evolution 
of Ca, we can further define £„(/) as the unique zero-average solution of Eq. (44), considered 
for fixed values of Ca, i.e. 

/f dl 
dlGa{l]Ca)\c^=c4i) = —Fa{l;ca). (45) 
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The indefinite integral in Eq. (45) is defined as the unique zero-average periodic primitive 
of the zero- average (periodic) function Ga{l)- During that integration, the arguments Ca are 
kept fixed, and, after the integration, they are replaced by the slowly drifting solution of 
the averaged system, given by Eqs. (43). Note that Eq. (44) yields Cq, = 0{Ga), which was 
assumed above, thereby verifying the consistency of the (approximate) two-scale integration 
method used here. 

We are now in a position to apply the above described method of variation of arbitrary 
constants, which gave us the evolution equations for Cq, and Cq,, to GW phasing. We use 2PN 
accurate expressions for the dynamical variables r, f, and entering the expressions for 
/ix and given by Eqs. (6). To do the phasing, we will solve the evolution equations for 
{ci, C2, Q, Ca}, given by Eqs. (39), on the 2PN accurate orbital dynamics, given in Eqs. (26). 
This leads to an evolution system, given by Eqs. (43) and (44), in which the RHS contains 
terms of order 0{c"^) x [1 + 0{c-^) + 0(c"^)] = 0{c~^) + 0{c~^) + 0{c~^). In the next 
section, as a first step, we will restrict our attention to the leading order contributions to 
Ga and Ga, which define the evolution of {ca, Cq,} under gravitational radiation reaction to 
0{c~^) order. We then impose these variations, via Eqs. (32) and (33), on to /ix and 
given by Eqs. (6). This will allow us to obtain gravitational wave polarizations, which are 
Newtonian accurate in their amplitudes and 2.5PN accurate in orbital dynamics. We name 
the above procedure 2.5PN accurate phasing of gravitational waves. Since G^'s create only 
periodic 2. 5PN corrections to the dynamics, in this paper, we will not explore its higher PN 
corrections. However, in a later section, we will present the consequences of considering PN 
corrections to Ga by computing (9(c^^) contributions to relevant This is required as 
Ga directly contribute to the highly important adiabatic evolution of and 

Up to now we have assumed, for concreteness, that the two constants ci and C2 were the 
energy and the angular momentum, respectively. However, any functions of these conserved 
quantities can do as well. In view of our use of the generalized quasi-Keplerian representation 
to describe the orbital dynamics, it is convenient to follow [31] and to use as Ci the mean 
motion n, and as C2 the time-eccentricity et. This can be done by employing the 2PN 
accurate expressions for n and Cj in terms of £ and J (or rather E and j), derived in [29, 30]. 
Firstly, this will require us to express 2PN accurate orbital dynamics in terms of /, n and 
Secondly, using n and Ct, instead of £ and J , as ci and C2, we need to derive the evolution 
equations for and ^ in terms of /,n and Cf. This will follow straightforwardly 
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from Eqs. (35). Using these expressions, the evolution equations, namely Eqs. (43) and 
(44), for {n, et, ci, c\, h, ej, ci, c\} will be obtained in terms of /, n and et. 

As mentioned earlier, we restrict in this paper the conservative dynamics to the 2PN 
order. Below, we present the 2PN accurate orbital dynamics, given by Eqs. (32), explicitly 
in terms of {l,n,et). This straightforward computation employs explicit expressions for the 
orbital elements of generalized quasi-Keplerian representation, in terms of E and j available 
in [29, 30]. The relations we need are: 



ar{n,et) 



Gm 



n 



1/3 c f2/3 f4/3 

l-^(9-r/) ■ ^ 



(72 + 75 r/ + 8 r/^ 



72 

^-1^(360 - 144 ry) - J^^i^^^ " 198 r/) 



(46a) 



^2/3 

er{n,et) = et{l + —— (8 - 3 77) + 



^4/3 



288 - 242 7/ + 21 r/2)ei2 



2 ' " ■ 24(l-e2)3/2 

390 - 308 r/ + 21 r/^^ ^1 - + (180 - 72 r/) (1 - e 

^4/3 



(46b) 



n,et) = et{l + e^^i-v) 



1152 - 656 r/ + 41 7/2)6*2 



96 (1 - e2)3/2 

1968 - 1088 7/ - 4 r/2^ ^1 - el + (720 - 288 r/) (1 - e: 



k{n,et) 



^4/3^4 



+ 



^4/3 



8v^ 



4(l-e2)2 



:(4 + 7/)7/et 



(51 - 26 7/) 6^ + (78 - 28 77) 



K5-277), 



2^^(1-377)7/6?, 
^e''c' 2 3 

32 (1 - 6,2)2 ' 



(46c) 
(46d) 
(46e) 
(46f) 

(46g) 
(46h) 



where ^ = 



Gmn 



We note that the generalized quasi-Keplerian orbital elements, given in 



terms of E and j in [29, 30, 49], can easily be expressed in n and et using following 2PN 
accurate relations for —2E and —2Ej^, which read 



~2E = e^^ii 



^2/3 


15 — r/ 


+ 


^4/3 


12 




24 


M20 48 


77^ 


1 


t \ 
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15 — 15 ?7 — r/" 



(47a) 



-2Ef = (1-e 



^2/3 



4(1 



(17 -7r/) et^ + 9 + r] 



'A? 



(360 - 144 7]) et^ \/\-e\ + (225 - 277 ^ + 29 if) et^ 



24(1 -e 

■ (210 - 190 r/ + 30 rf) tt^ + 189 - 45 r/ + r/' 



(47b) 



These two relations easily follow from inverting the 2PN accurate relations for the orbital 
period P = ^ and e\ in terms of E and j presented in Eqs. (19) above [ See [29, 30] ]. 
In addition, to compute expressions for r and 0, we use the following relations 

du 



dS_ 

'm 

dW 

"m 

du 
'dl 

dv 
du 



ttr Cr sm u 



dl 



1 - ej cos n 1 + ^ (/t cos v + gt) 



du 



(1 - eiy/' 



(48b) 
(48c) 
(48d) 



1 — cos n 

The radial motion, defined by r(/, n, Ct) and r(/, n, ej, reads (both in the compact form and 
in 2PN-expanded form) 



r = S {I, n, et) = ar{n,et) {1 — er{n,et) cos u) 



(^)V3(i-e, cos«)a- 



^2/3 



(6 — 7 ?7) cos m + 18 — 2 77 



^4/3 



72 v/(r^^ef)3(l -et cosn 
-(72 + 75 7/ + 8 7/2)6*2 - 234 + 273 r/ + 8 r/ 

-36 (1 - e^) (5 - 2 r/) (2 + cos n) 



6 (1 — e* cosn) [ 

(72 - 231 7/ + 35 7/2) (1 - 6*2) e* cos u 



(49a) 



dS,, , (GmnY 

n—{l,n,et) = 7- re* smw 

dl (1 — e* cosm) 



(6-77/) 



^ (- (72 - 231 7/ + 35 7/2) (e* cosm) 

cosn)'^ L V 



^4/3 

(1-e* 

+ (216 - 693 7/ + 105 r/2) (e* cos uf + (324 + 513 r/ - 96 7/2) e* cos n 
-(36 + 97/) 7/6*2 - 468 - 15 r/ + 35 7^2 
36 / 

(1 — e* cos m) (4 — e* cos n) (5 — 2 r/) 



(49b) 
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In the above equation, the eccentric anomaly u = U{l,n,et) is given by inverting the 2PN 
accurate Kepler equation, Eq. (30), connecting / and u, i.e. in explicit form 

^4/3 



u — et sm u 



I et sin u \/\- e\ i] (4 + 77) 

t cosn) [ 



+12 (5 - 2 77) (n - v) (1 - et cosn) 
The angular motion, described in terms of </> and </>, is given by 



(50) 



0(A,O = A + 1^(0 



(Z) = V — uAr et sin u + 



(51a) 



(1 - e?: 



— n + sin u 



^4/3 



I 4 Vl - e? {\-et cosm)' (^|- 



(102 



32(l-e?)5/2 (1 _e, cosn)3 
-52 7/)et2- 156 + 56 7/|ej cosu + r/ (4 + r/) + (l02 - 60 - 2 r/^) e*^ 

+ 156 - 52 r/ + r/2^ + (l - e^^) (^((3 e^^ + 12) r/ - 8) (e* cosn)^ 

+ ((8 - 6 r/) et^ + 8 - 24 r/) {et cos m) - 12 r/ e*^ - (8 - 27 r/) Ct^^ 



Ci sin n + (1 — cos u) 



48 (1 - e\f (5 - 2 r/) - 8 ( (51 - 26 r])et 



-78 - 28 7/ 



(1 — et cos n) 



(51 - 26r/)et^ 



+78 



28 7/^ Vl - - (30 - 12 r/) (1 - 2 e^^ + e^) w| 



(51b) 



n 



^2/3 



(1 — et cosM 
- (4 - r/) + 3 



'^^^^ (1 _ gS-) (^]^ _ cosm) 



i^l^ 1 

12 (1 — cosm) 



(1 — 77) et cosM 
1 



(1 - e2)3/2 



18 (1 — cos u) 



[et cosu - 26^2 + 1) (5 - 27/)| + (-(9 " 19^ - 14r/^ 

-36 + 2r/-8r/2 ) (e* cosn)^+ ( -(48 - 14 r/ + 17r/^)e/ 



(69 - 79 r/ + 4 r/2) Cf^ + 114 + 2 r/ - 5 r/^^ (e* cosn) 
(^-(6 - 32 7/ - 7/2)e/ + (93 - 19 7/ + 16 t/^) e^^ 
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-222 + 50 r/ + j (g^ cosm) - 6r/ (1 - 2r/) 

+ (54 - 28 r/ - 20 r/^) e*^ - (l53 - 61 r/ - 2 r/^) e*^ + 144 - 48 r/ 



}] } ■ ^^^^ 



The explicit form of above follows from Eq. (32b). 



In addition to the above explicit expressions, we also need to evaluate the RHS of 
Eqs. (35a)-(35b) and, in particular, the 2PN accurate partial derivatives of n and with 
respect to the relative velocity v. To get these, one could combine Eqs. (19) with the expres- 
sions for E and j in terms of relative position and velocity, rather than in terms of position 
and momenta as is usual in the ADM formalism. To the desired 2PN order, one may either 
start from the ordinary Lagrangian L(x, v) in ADM coordinates ( See [45] for the explicit 
construction of this Lagrangian ) or (simply) by inverting the basic Hamiltonian equation 
V = dH/dp, to get v in terms of p. However, in the next section, we require expressions for 
E and j only to the well known Newtonian order. 

V. 2.5PN ACCURATE PHASING 

Let us recall that our method is general and can be applied, in principle, to any PN 
accuracy. For instance, we could study the effect of the 0{c~^) + 0{c~'^) radiation reaction 
on the 3PN conservative motion. However, in this work, we limit ourselves, for simplicity, to 
considering the effect of the 0{c~^) radiation reaction on the 2PN motion. Accordingly, we 
shall, each time it is possible, truncate away all effects that would correspond to the (9(c~^) 
level or beyond. As we shall see, this approximation is probably sufficient for oscillatory 
effects (in the sense of the decomposition, given in Eq. (40) ), which are the primary 
focus of this paper. We shall discuss below, how our method also justifies the usual way 
of deriving the secular effects linked to the radiation reaction, and we shall obtain more 
accurate expressions for them. 

This section begins by providing inputs necessary for computing evolution equations for 
the set {cq,Cq,}, to the 2.5PN order, where the index a = n,et,ci, and c^. As just said, 
we require A' to 2.5PN order for this purpose. The 2.5PN expressions for A' will have to 
be in the ADM gauge, as our conservative 2PN dynamics is given in the same gauge. The 
expression for relative reactive acceleration, to 2.5PN order, in the ADM gauge, available in 
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[51], reads 



where v"^ 
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I2v' - 15 r' 



Gm 
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(53) 



r . At this point a nice technical simplification occurs. Though our 



formahsm consistently combines a 2PN-accurate, precessing motion with 2.5PN radiation 
reaction, the RHS's of Eqs. (39) are technically given by the product of 0{c~^) reaction 
terms by orbital expressions given as explicit PN expansions 0(c°) + (9(c~^) + (9(c~^). 
Therefore, if we decide, in a first approach, to neglect 0{c~^) contributions to the phasing, 
we can simplify the RHS's of Eqs. (39) by keeping only the leading terms in the orbital 
expressions. This formally means that it is enough to use Newtonian-like approximations 
for all orbital expressions appearing in Eqs. (39). For instance, we can simply use r ~ 
(GM/n2)V3(i_e^cosM), n ~ {-2Ei)sfl'^/Gm^ [{2Gm/r - v'^)f''^ /Gm, etc. in Eqs. (39). 
Note, however, that this does not at all mean that we are approximating the orbital motion 
as being a non-precessing Newtonian ellipse. In all expressions where they are needed, we 
must retain the full PN expansion. For instance, in the contribution A, given by Eq. (63b), 
(see below) to = \ + W{l), we must keep the 2PN accuracy for the precession rate n{l + k), 
and augment it by the effect of the time variation of n{t) and k{n(t), et(t)), as discussed there. 

Finally, the leading evolution equations for {^, ^} in terms of u{l, n, et),n, and 

Ct, follow as 
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where x = (1 — e* cosm), ^ = and n = u{l,n,et). We are in a position to explore 

the secular and periodic variations of {cq} to 0{c~^), which will be done in the next two 
subsections. 



A. Secular variations 



Using the above set of equations with Eqs. (43) and (44), we obtain the differential 
equations for {cq,,Cq,}, where the index a = n,et,ci,cx- Let us first consider the secular 
variations of Cq, given by Eqs. (43). One remark is that, after using an /-variable formulation 
to separate the secular variations from the oscillatory ones, we can, at the end, re-express 
the secular result, Eqs. (43), in terms of the original time variable t. This leads to 

^ = F^ic.) (55) 

where is the /-average of the RHS of the t- variation of the c^'s, see Eq. (34): Fa{ca) = 
(27r)~^ J^^ dl Fa{l; Ca)- Among the secular variations, let us first discuss the secular variation 
of Q and c\. The 'source term' for the secular variations of q and c\ is the /-average of Fi 
or Fa, i.e. modulo a (secular) factor n{ca), the /-averages of Gi = Fi/n, or G\ = Fx/n, i.e. 
the /-average of the RHS's of Eqs. (54c) and (54d). A look at the RHS's show that, being 
of the form sinu /(cosu), they are odd under u —u, so that their average over d/ ^ (1 — 
etCOSu)du exactly vanishes: Gi = = G\. In fact, this remarkable finding follows from the 
time-odd character of the perturbing force A', and therefore would also hold if we considered 
the radiation reaction to the accuracy 0{c~^) + 0{c~'^) + (9(c~®) + 0{c~^). [ We stop at 
order because of the conceptual subtleties arising in the meaning of radiation reaction at c~^° 
order, which is the first order where nonlinear effects linked to the leading 0{c~^) radiation 
reaction enter]. Note that the 0{c~^) contribution to radiation reaction comes from the 
tail contributions, which in its exact form is given by an integral over the past [52]. This 
correction is time-reversal asymmetric without being simply time-reversal antisymmetric. 
However, when approximating that integral as a function of the instantaneous state, it 
becomes a time-reversal antisymmetric function of x and v [53]. Indeed, Ci and C2 being 
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even under time-reversal, the partial derivatives dci/dv^, dc^jdv'^ ^ in Eqs. (35) are time-odd. 
When is time-odd, Eqs. (35) then imply that dci/dt and dc2/dt are time-even. Then in 
Eq. (35c), dS/dl is time-odd (because r is even, but / is odd), dS/dca is even and dca/dt is 
time-even, so that dci/dt ends up being time-odd. The same conclusion is found to hold for 
dcx/dt, thereby ensuring the absence of secular variations for both q and cx- 

^ = ; ci{t) = ci{0), (56a) 
^ = ; c,(t) = c,(0). (56b) 

Turning now to the secular variations of n and Ct, Eq. (55), we note that they reduce to 
the usual 'adiabatic' estimate of the secular variation of constants, namely 

dca ^ <9ca(x,v) ,i 
^ = ^^^^ ^' 

where ()/ denotes an average over an (instantaneous) orbital period. If we were using as 
Ci and C2 the system's dimensionless energy E and angular momentum j, Eq. (57) is the 
usual way of estimating the secular change of E and j under the influence of a perturbing 
acceleration A'. Applying Eq. (57) to the case where Ci = n and C2 = et is easily seen to lead 
simply to a coupled differential system for n{t), et{t), which is strictly equivalent (under the 
map n = n{E,j), et = et{E,j)) to the just mentioned secular evolution system for E and j. 
The /-average of the RHS of Eqs. (54a) and (54b) in the leading 0{c~^) approximation, as 
mentioned earlier, only lead to the leading terms in the secular evolution of n and e^. Since 
the RHS's of Eqs. (54a) and (54b) are expressed in terms of u, it is convenient to do the 
orbital average by expressing it as an integral over using u, using (i/ ^ (1 — et cosu)du. The 
resulting definite integrals may be easily computed, using [54], which give 

_1_ f^^ du _ 1 / 1 \ 

211 J, (1-e, cos«)^+i " (l-e?)(^+i)/2 ^^^/T^J' ^^^> 

where P/v is the Legendre polynomial. Using Eq. (58) in Eqs. (54a) and (54b), we obtain 
leading 0{c~^) corrections to ^ and From these expressions, using dl = ndt, we obtain 
^ and which read 

at dt ' 

dn _ (G-)^/^-"/^^/96 + 292e,^ + 37e,4+0(l), (59a) 
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These results are equivalent to the old results of Peters [4] based on balance argument 
between far-zone fluxes and local radiation damping. In the next section, we will present 
2PN accurate expressions, providing corrections to (9(c~^°) for ^ and 

Let us finally note that formally one can analytically solve the coupled evolution system 
by successive approximations, reducing it to simple quadratures. For instance, at the leading 
order where one keeps only the 0{c~^) contributions, one can first eliminate t by dividing 
dn/dt by det/dt, thereby obtaining an equation of the form dlnn = fo{et)det. Integration 
of this equation yields 

where is the value of et when n = rij a result first obtained by Peters in [4]. 

Inserting Eq. (60) into the leading evolution equation for et then leads to an evolution 
of the form det/dt = goi<^t) which can be done by quadrature: t = J detg^^iet) + constant. 
We can then insert back the leading result, Eq.(60), into the leading correction terms of 
the evolution equation for dlnn/det to get again a decoupled equation of the form dlnn = 
f2{et)det which can be integrated. Continuing in this way would give the function n{et) in the 
form of a expansion, which would lead to an explicit decoupled equation for the temporal 
evolution of et: det/dt = g{et) = go + c~'^g2, which can again be solved by quadrature. 
This procedure may easily be extended to 0{c~^) order. At the leading 2.5PN order, we 
have checked that the temporal evolution for (n, ej, obtained by solving coupled differential 
equations, Eq. (59), is in excellent agreement with those given by the above mentioned 
procedure. 



B. Periodic variations 



Let us turn to the differential equations, which give at 0{c ^) order, orbital period 
oscillations to our dynamical variables. They read 
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where n and 6*, on the right hand side of these equations, again stand for n and e*. Here the 
RHS's of Eqs. (61) are zero-average oscillatory functions of /. [The RHS's for Eqs. (61c) and 
(61d) are in fact identical to the ones of Eqs. (54c) and (54d), in view of our previous result 
= Ga = 0, except for the fact that they are expressions in terms of n and e*, instead of 
n and 6*.] 

One can analytically integrate Eqs. (61) to get n, e*, c/, c\ as zero average oscillatory 
functions of /. We find, when expressed in terms of u 
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where Pt = — — -■ Finally, let us consider the way the previous results feed in to the basic 
angles / and A entering our perturbed solution Eq. (24). From the definitions for l{t) and 
A(t), as given in Eqs. (33), we see that we can also split these angles in 'secular' pieces, say 
I, A and 'oscillatory' ones /, A, as 



m = T{t) +l{l-Ca{t)) 

\{t) = \{t) + A(/; Cait)) , where 
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We note, based on earlier results, that ci{t) = 6; (to) and c\(t) = Cxito) are constants. The 
'oscillatory' contributions to / and A are given by 
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Here k = idk/dn)n + ((9/c/(96*)e* denotes the oscillatory piece in k and / dlfil) denotes the 
unique zero average primitive of the zero-average periodic function of /, /(/). 

To complete our study of the Oic~^) oscillatory contributions to the phasing, we see from 
Eq. (64a) that we need to integrate n/n and add it to the previous result for q(/). [Note 
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that they appear together in the phasing formula.] We find 

mca) = Y^^^^^^|(602 + 673et')x+ (314- 203et'-llle3lnx 
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where Ci{l) is given by Eq. (62c) and pt = ^^—3. The exphcit expression for A at the 
2.5PN order is simply given by Eq.(65) with q(/) replaced by the expression for with cx{l), 
given by Eq. (62d). This is so since contributions to A(/) arising from the periastron advance 
constant k appears at 0{c~'^). In the next subsection, we plot analytic results obtained in 
these subsections and their infiuences on hx and h+. 



C. Graphical representation of the results 

We begin the subsection by illustrating the temporal evolution of Ca and c^. Next, we 
show the combined effects of these secular and periodic variations on basic angular variables 
that appear in the expressions for and h^. Finally, we exhibit and h+ evolving 
under gravitational radiation reaction and point out various features associated with post- 
Newtonian accurate orbital motion. In these figures, we terminate the orbital evolution 
when j2 = 48. This criterion, as explained earlier, is chosen to make sure that the orbit 
under investigation is a shrinking slowly precessing ellipse. We have also taken advantage 
of the 'scaling' nature of the problem to plot only dimensionless quantities in terms of 
dimensionless variables. The conversion to familiar quantities like orbital frequency / (in 
Hertz) is given by / = ^ = 2^ = 3.2312 x 10^.^ This implies that for a compact 
binary with m = Mq and ^ = 10"^^ the orbital frequency will be ~ 30 Hertz. 

In Figs. (2) and (3), we plot n/rii (where rii is the initial value of n), h/n,et,et,ci,ci, c\ 
and Ca as functions of which gives evolution in terms of elapsed orbital cycles. We clearly 
see an adiabatic increase (decrease) of n (e*) as well as the periodic variations of n and e*. 
As expected, we also observe no secular evolution for q and cx, but clearly see periodic 
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variations in c/ and cx. In order to illustrate the effect of the secular and periodic variations 
in the above orbital variables on the basic angles / and A, in Fig. (4), we plot scaled dl/dt 
and d{X — l)/dt as functions of Thefi gure shows periodic oscillations superposed on the 
slow secular drift. Finally, in Figs. (5) and (6), we plot scaled /i+ and hx as functions of 
We employ, for these figures, polarization amplitudes, which are Newtonian accurate while 
the orbital motion is 2.5PN accurate. We clearly see 'chirping' due to radiation damping, 
amplitude modulation due to periastron precession and also orbital period variations. 

Using the scaling argument mentioned earlier, we note that Figs. (2)- (6) may be used 
to illustrate the various aspects of a compact binary inspiral from sources relevant for both 
LIGO and LISA. For instance, if we choose m = 2.8Mq, the variation of ^ from 2.069 x 10"^ 
to 3.0186 X 10~^ in around 227 orbital cycles corresponds to orbital frequency variation from 
150 Hz to ~ 217 Hz in ~ 8.15 seconds. Similarly, the choice m = IO^Mq corresponds to 
a binary inspiral involving two supermassive black holes, where orbital frequency increases 
from ~ 4.2 x 10"^ Hz to ~ 6.2 x 10"^ Hz in ~ 2.7 days. 

Finally, we note that Figs. (2)- (6), were drawn mainly to exhibit more clearly the exis- 
tence of periodic variations in orbital elements (analytically investigated for the first time in 
the present paper) due to the radiation reaction. Evidently, as these variations scale as 
see Eqs.(62), they become quite small if the binary is not near the plunge boundary. Our 
work is important, even when the binary is not near the LSO, as it shows, in a technically 
clear manner, how to describe the exact phasing as the sum of the usually considered 'adia- 
batic' phasing (involving only secular variations) and a normally neglected 'post-adiabatic' 
phasing (involving only the relatively smaller 'periodic' variations). More about this in the 
conclusions below. 

VI. PN ACCURATE ADIABATIC EVOLUTION FOR n AND et 

One of the useful results of the present work are Eqs. (56) and (57). Indeed, these 
equations provide a clear justification for the usually considered 'adiabatic' approximation 
(in cases, where one is sufficiently away from the plunge boundary so that one can safely 
neglect the additional periodic contributions and treat orbits to be quasi-Keplerian) . More 
precisely, we have, earlier, proved that Eqs. (56) would be valid, even if we are considering 
radiation reaction to the accuracy 0{c~^) + 0{c~'^) + 0{c~^) + 0{c~^). Concerning Eqs. 
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(57), following its derivation again, we see that (if we were to define Ca to sufficient accuracy, 
by considering the conserved quantities of the 'conservative part' of the dynamics) it is valid 
up to terms which are quadratic in the radiation reaction, i.e up to 0{c~^^). Therefore, to 
those very high accuracies, our work shows that the secular part of the phasing (in the sense 
of the decomposition, as in Eq. (40) ) is essentially given by the simple averaged result, 
given by Eqs. (56) and (57). Then, as usual, we expect that the averaged losses of the 
mechanical eneigj and angular momentum of the system, appearing in Eqs. (57), are equal 
to the corresponding far zone (FZ) fluxes of energy and angular momentum in the form of 
radiated gravitational waves. 

To complete this work, let us briefly show how to obtain, in ADM coordinates, the 2PN 
accurate secular changes in n and Cj, equivalent to corresponding 2PN accurate FZ fluxes 
of energy and angular momentum. The differential equations for n and et are computed, 
following the PN accurate calculations presented in [17]. These computations require PN 
corrections to orbital averaged expressions for the far-zone energy and angular momentum 
fluxes and PN accurate expressions for n and ej, all of these expressed in terms of orbital 
energy and angular momentum. The PN accurate expressions for ^ and ^ are obtained 
by differentiating PN accurate expressions for n and Cj w.r.t. time and then heuristically 
equating the time derivatives of orbital energy and angular momentum to orbital averaged 
expressions for the far-zone energy and angular momentum fluxes. For the ease of imple- 
mentation, we split 2PN accurate computations of ^ and ^ into two parts. The first part 
deals with the purely 'instantaneous' 2PN corrections and the second part considers the so 
called 'tail' contributions, appearing at the 1.5PN (reactive) order [55]. The computations to 
get 'instantaneous' contributions begin with 2PN corrections to far- zone fluxes, in harmonic 
gauge, in terms of r, r and available in [28]. Using 2PN accurate relations connecting the 
dynamical variables in harmonic and ADM coordinates, given by Eqs. (Al) in Appendix A, 
we obtain expressions for the far- zone fluxes in ADM coordinates. These far-zone fluxes are 
orbital averaged, using 2PN accurate generalized quasi-Keplerian parametrization for ellipti- 
cal orbits, following lower PN computations done in [17]. We then compute time derivative 
of PN accurate expressions for n and Ct and equate resulting time derivatives of orbital 
energy and angular momentum to orbital averaged expressions for the far-zone energy and 
angular momentum fluxes respectively, to get PN accurate ^ and ^ in terms of E,j,m 
and 7]. Finally, we use Eqs. (47) to obtain the differential equations for n and et in terms of 
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n, et, m and r]. The tail contributions to ^ and ^ are already available, in slightly differ- 
ent forms, in [17] and we only rewrite these expressions in our n and ej variables. Adding 
these 'instantaneous' and 'tail' contributions gives us the following 2PN accurate evolution 
equations for n and et- 

dn [^n , ^ipn , ^i.spn , ^2Pn1 v 

* = gmV^ 1" /■ 

'Bl . _(^!^^gl££|e» + e,'™ + e,'-^™+e,^™,} (66b) 

where n^, fi^™n'™ , Cf^, e/^'^ and e*^™ denote 'instantaneous' contributions to 2PN order, 
while n"^^^ and Cf^'^^^ stand for the 'tail' contributions. Using 2PN accurate orbital 
representation and far-zone energy and angular momentum fluxes, we have computed 2PN 
accurate 'instantaneous' contributions to ^ and whose explicit forms are given by 
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where, as in earlier instances, n and et on the right hand side of these equations stand for n 
and et. 

The 'tail' contributions to ^ and which appear at 1.5PN order, are derivable using 
Keplerian orbital parameterization and 'tail' corrections to orbital averaged expressions for 
the far-zone energy and angular momentum fluxes available in [53, 56]. For the ease of 



presentation, we display below 'tail' contributions to ^ and rather than ^ and ^ 
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where both ke and kj are expressible in terms of infinite sums involving quadratic products of 
Bessel functions Jpipet) and its derivative J'^ipet). For completeness, the explicit expressions 
for ke and kj, given in [53, 56], are listed below 



p=i 



1 1 1 2/1 3 



2 I 3 



4 7 

-- + - 

et 



3et 



Jp{pet)J'^{pet) + {J'p{pet)f 



2 + el 



+00 2 



(69a) 





2 











Jpipet) J'pipet) + 2p(^l - (j;(pet))2| . (69b) 



We have checked that to IPN order the above equations are consistent with expressions for 
%^ and computed in [17]. At 2PN order, the above expressions are also consistent with 



corrected formulae for ^ and available in [28, 32] 
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VII. CONCLUSION 



Let us summarize what we have proposed in this paper and point out possible extensions. 
We have provided a method for analytically constructing high-accuracy templates for the 
gravitational wave signals emitted by compact binaries when they move in inspiralling slowly 
precessing eccentric orbits. In contrast to the simpler problem of modeling the gravitational 
wave signals emitted by inspiralling circular orbits, which contain only two different time 
scales (orbital period and radiation reaction time scale), the case of inspiralling eccentric 
orbits involves three different time scales: orbital period, periastron precession and radiation- 
reaction time scales. An improved 'method of variation of constants' is used to combine 
these three time scales, without making the usual approximation of treating adiabatically 
the radiative time scale. By going to a suitable center-of-mass frame, the transverse traceless 
(TT) radiation field and hence the GW polarizations are expressed as PN expansions of the 
form given in Eq. (4) in harmonic coordinates involving only the relative position and 
velocity. The polarisations can be rewritten in terms of the ADM positions and velocities by 
using the contact transformations available to move from the harmonic coordinates to the 
ADM coordinates. In the ADM coordinates, the unperturbed 2PN (3PN) motion may be 
explicitly solved by a generalized quasi-Keplerian representation involving two angles / and 
A and four constants of the 2PN (3PN) motion Cq = Ci, C2, Q and c^. By a Lagrange method 
of variation of constants, this unperturbed solution is used to prescribe a general solution 
to the perturbed (by radiation-reaction) system of the form given by Eqs. (32) and (33), in 
terms of the varying c^'s. Among the four c^'s, two (q and cx) are found to be constants, 
while the other two c^'s satisfy two coupled first order differential equations. A two scale 
decomposition of the c^s is made to model the combination of the slow (radiation reaction, 
secular) drift and the fast (orbital time scale, periodic) oscillations. This allows us to 
decompose the TT gravitational wave amplitudes or polarizations into a part associated with 
the secular variations and another part associated with the fast oscillations. If one expands 
in the fast oscillations, the oscillatory contributions to the phasing would be equivalent to 
new additional 0{v^/c^) contributions (which stay always small ) to be added to the usually 
considered 2.5PN amplitude hfj in the adiabatic approximation. Therefore, as long as an 
amplitude correction 0{v^/c^) is not needed our results show that it is enough to use only 
secularly varying n and e^. Note that our description of secular effects automatically include 
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an effect of a secular acceleration of the periastron precession, analogous to the usual secular 
acceleration of the orbital motion. Namely, the secular part of the angle A — / measuring 
the periastron longitude varies as 



A — / = / kndt 



(70) 



where k{t) = 5$/27r is the secularly varying periastron precession per orbital period. 
Schematically, thus our work may be summarized as follows: 



h, /^harm J,harm -harm J,harm 



'j+.x 0*™, r*™'', 0*"'') 




h+,x{n, et, l,X} + ^ f^/i^.xl". et, l A) 



In this paper, the orbit of the binary was treated to be an inspiralling, slowly precessing 
ellipse, which prevented us from approaching the LSO. However, as mentioned earlier, using 
the EOB approach, we intend to explore the orbital dynamics and hence the evolution of 
gravitational wave polarizations near the LSO in the near future. There are also quite a few 
generalizations which can be tackled using the formalism presented here and we list a few 
of them here. In this paper, the conservative dynamics was restricted to the 2PN order and 
therefore, a natural extension will be to incorporate the 3PN conservative dynamics. We 
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also restricted our approach to compact binaries consisting of non-spinning point masses. It 
is possible to extend the generalized quasi-Keplerian parametrization and hence our method 
to spinning compact binaries. To do that, first, one needs to extend the quasi-Keplerian 
representation to include the effects due to spin-orbit and spin-spin interactions, which 
requires generalizing a restricted analysis done in [57] and this is under investigation [58]. 
Finally, starting from the gravitational wave polarizations from inspiralling eccentric binaries 
in the time-domain, it should be possible to perform a spectral analysis and see how their 
power spectrum depends on various orbital elements like n, Ct and i. 
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APPENDIX A: THE CONSTRUCTION OF PN CORRECTIONS TO /ix AND h+ 

In this appendix, we sketch the procedure to compute PN corrections to and in 
ADM coordinates, in terms of the dynamical variables r, r, and 0. It is clear from Eqs. 
(1), (3) and (6) that PN corrections and /ix require PN corrections to hj^ . The 'in- 
stantaneous' 2PN accurate contributions to hj^ in harmonic coordinates in terms of the 
components of n and v, r, r and v'^ were computed in [28]. To obtain similar expressions 
in ADM coordinates, we employ the 2PN accurate contact transformations linking the har- 
monic and ADM coordinates, given in [45] , which prescribe the way to relate the dynamical 
variables in these coordinates. We list below the transformation equations relating x and v 
in the harmonic coordinates to the corresponding ones in ADM [32]: 



Gm 

Xh = xa + 



CjfTfl 

(5i;2 - r^) r/ + 2 (1 + 12?;) 



-18?7rr v} , (Ala) 
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77 + 4 
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13 1;' 



17f2 - 42 



2 U. 

r J 



(Alb) 



The subscripts 'H' and 'A' denote quantities in the harmonic and in the ADM coordinates 
respectively. Note that in all the above equations the differences between the two gauges 
are at the 2PN order and hence no suffix is used for the 2PN terms. We do not list the 
transformation relations for r, f and as they easily follow from Eqs. (Al), as r = |x|, r r = 



X ■ V and v 



V ■ V 



^2 _|_ ^2 02_ 



Using Eqs.(Al), the 2PN corrections to hjj in ADM coordinates can easily be obtained 



from Eqs. (5.3) and (5.4) of [28]. For economy of presentation, we write {hf )a in the 



following manner, {hj^ )a = {hfj 



Gorrections\ where (/i^ )a represent the metric 



perturbations in the ADM coordinates. {hJj)o is a short hand notation for expressions 
on the R.H.S of Eqs. (5.3) and (5.4) of [28], where N, n, v, w^, f, r are the ADM variables 
Na, Ha, Va, t'A5 ''^A, ''"A respectively. The ^Gorrections^ represent the differences at the 2PN 
order, that arise due to the change of the coordinate system, given by Eqs. (Al). As the 
two coordinates are different only at the 2PN order, the ^Gorrections^ come only from the 
leading Newtonian terms in Eqs. (5.3) and (5.4) of [28]. 
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We now have all the inputs required to compute the 2PN corrections to /i+ and /ix in 
ADM coordinates, in terms of r, 0, f and (p. We write, schematically, the expression for 



(hf^A as 



{h^j )a = a™ Vij + ann riij + Unv n{iVj) . (A3) 

We apply exactly the same procedure which gave us the Newtonian expressions to and 
hx from Newtonian contributions to /i^^. Since the explicit expressions for the final 2PN 
accurate 'instantaneous' GW polarization states are too lengthy to be listed here, we write 
schematically 



a,; 



+ (■••)cos2</) + (■■■)sin2</) 



+ 



cos 2(f) 
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sin 2 
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) sin 20 
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In the above and what follows (■ ■ ■ ) denotes various coefficients expressed in terms of 
r, 'r,0, mi,m2 and i. The structure of the PN expansion of the coefficients aij above is 
the following: 



a„, ~ 1 



(■■■)cos0 + (■■■) 



sm I 



+ 



(■ ■ ■ ) cos 30 + (■ ■ ■ ) sin 30 



(■■■) + (■ • ■) COS20+ (■ ■ ■) sin20 
(■■■) + (■■■) cos 40 +(■■■) sin 40 



+ 



(A5) 



ann and have similar expansions with the exception that for anv the leading order term 
is at ^ order. We note that, similar to Eq. (6), the above sketched expressions for /i+ and 
/ix are in a form suitable to apply our phasing formalism. 
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FIG. 1: The 2PN accurate effective radial potential Wj{u) as a function of the dimensionless radial 
variable ^ = for various values of the dimensionless angular momentum j. The point, marked 
c, denotes a stable circular orbit, while the line, noted e, stands for a precessing elliptical orbit. 
The line with label p denotes an elliptical orbit which is about to plunge. Note that the left end of 
the line p is tangent to the effective potential, and corresponds to an unstable circular orbit. The 
plots are for t] = 0.25. 
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These variations are governed by the reactive 2.5PN equations of motion. Periodic nature of the 
variations in n and Cf are clearly visible, and e{ denote initial and final values for the time 
eccentricity et, while and ^-^ stand for similar values of the adimensional mean motion '^^^ . 
The plots are for rj = 0.25 and the evolution is terminated when = 48. 
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FIG. 3: The plots ci,ci,cx and ca against orbital cycles, given by Similar to Fig. 2, these 
variations are governed by the reactive 2.5PN equations of motion. Periodic nature of the variations 
in ci and c\ as well as the constancy of q and cx are clearly visible. The symbols have the same 
meaning as in Fig 2. 
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FIG. 4: The plots showing scaled time derivative of I and A — / as function of orbital cycles, given 
by The right panels show clearly both the secular drift and the periodic oscillations of the 
plotted quantities. Similar to earlier figures, these variations are governed by the reactive 2.5PN 
equations of motion. The initial and final values of the relevant orbital elements are marked on 
the plots. The plots are for r] = 0.25 and rii is the initial value of the mean motion n. 
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FIG. 5: The scaled h+ ( Newtonian in amplitude and 2.5PN in orbital motion) plotted against 
orbital cycles, given by The 'chirping', amplitude modulation due to periastron precession are 
clearly visible in the upper panel. In the bottom panel, we zoom into the initial stages of orbital 
evolution to show the effect of the periodic orbital motion and the periastron advance on the scaled 
The initial and final values of the relevant orbital elements are marked on the plots. The 
plots are for rj = 0.25 and the orbital inclination is i = ^. 
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FIG. 6: The plots for the scaled hy_ (Newtonian m amplitude and 2.5PN in orbital motion) as 
function of where I is the mean anomaly . The slow 'chirping', amplitude modulation due to 
periastron precession are clearly visible in the upper panel. In the bottom panel, we zoom into 
the initial stages of orbital evolution to show the effect of the periodic orbital motion and the 
periastron advance on the scaled /ix (t). The initial and final values of the relevant orbital elements 
are marked on the plots. The plots are for a binary consisting of equal masses, such that i] = 0.25 
and the orbital inclination is i = ^. The orbital evolution, as in earlier cases, is terminated when 
P = 48. 
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